function diff = model_2(param,aux)
   
    %% Solve model in the replacement region as there are no turnover costs 
    region = 'R';

    %% Calculate roots in the natural wastage region
    a=-1/2.*param.sigma^2; 
    b= -(param.mu-1/2.*param.sigma^2+(1-param.alpha).*(param.zeta));%param.s.*param.lambda+
    c=+param.r+param.zeta;%+param.s.*param.lambda
    param.gamma_N_1 = (-b+(b.^2-4.*a.*c).^(1/2))./2./a;
    param.gamma_N_2 = (-b-(b.^2-4.*a.*c).^(1/2))./2./a;
    param.gamma_R_1 = param.gamma_N_1;
    param.gamma_R_2 = param.gamma_N_2;
    clear a b c
    
    % Create bounds for G=m_u/m_l
    if param.K>0 
        G_l     = 1+(param.c+param.K).*(param.r+param.zeta)./(param.omega_0);%param.s.*param.lambda+
        G_u     = (1-1./param.gamma_N_1)./(1-1./param.gamma_N_2).*(1+(param.c+param.K).*(param.r+param.zeta)./(param.omega_0));%param.s.*param.lambda+    
        tmp_G   = fzero(@ (G) (param.omega_0 + rho_fun(0,region,param).*(param.c+param.K))./param.omega_0.*phi_fun(G,region,param) ...
                    -G.*phi_fun(G.^(-1),region,param), [G_l,G_u], optimset('display','off'));
        param.G = tmp_G;
        
        param.m_l = param.omega_0./(1-param.omega_1)./rho_fun(0,'R',param)./phi_fun(param.G,region,param);
        param.m_e = (param.omega_0+rho_fun(0,'R',param).*(param.c+param.K))./(1-param.omega_1)./rho_fun(0,'R',param)./phi_fun(param.G.^(-1),region,param);
        param.m_u = param.m_e;            
        
        param.J_1 = - (1-param.omega_1).*theta_fun(param.G,region,param).*param.m_l.^(1-param.gamma_N_1)./param.gamma_N_1./rho_fun(1,'R',param);
        param.J_2 = - (1-param.omega_1).*(1-theta_fun(param.G,region,param)).*param.m_l.^(1-param.gamma_N_2)./param.gamma_N_2./rho_fun(1,'R',param);
        param.m_h = fzero(@ (m) ((1-param.omega_1).*m./rho_fun(1,'R',param) - param.omega_0./rho_fun(0,'R',param) ...
            + param.J_1*m.^(param.gamma_N_1) + param.J_2*m.^(param.gamma_N_2) - param.c), [param.m_l,param.m_u], optimset('display','off'));
        
        if param.zeta==0
            param.B=(2.*(1-param.alpha)./param.sigma^2.*param.s.*param.lambda).^(-1)-log(param.m_h);
        end
        
    else 
        G_l=1+param.c.*(param.r+param.zeta)./(param.omega_0);%param.s.*param.lambda+
        G_u=(1-1./param.gamma_N_1)./(1-1./param.gamma_N_2).*(1+param.c.*(param.r+param.zeta)./(param.omega_0));%param.s.*param.lambda+    
        tmp_G= fzero(@ (G) (param.omega_0 + rho_fun(0,region,param).*param.c)./param.omega_0.*phi_fun(G,region,param) ...
                    -G.*phi_fun(G.^(-1),region,param), [G_l,G_u], optimset('display','off'));
        param.G=tmp_G;
        
        param.m_l = param.omega_0./(1-param.omega_1)./rho_fun(0,'R',param)./phi_fun(param.G,region,param);
        param.m_h = (param.omega_0+rho_fun(0,'R',param).*param.c)./(1-param.omega_1)./rho_fun(0,'R',param)./phi_fun(param.G.^(-1),region,param);
        param.m_u = param.m_h;

        param.J_1 = - (1-param.omega_1).*theta_fun(param.G,region,param).*param.m_l.^(1-param.gamma_N_1)./param.gamma_N_1./rho_fun(1,'R',param);
        param.J_2 = - (1-param.omega_1).*(1-theta_fun(param.G,region,param)).*param.m_l.^(1-param.gamma_N_2)./param.gamma_N_2./rho_fun(1,'R',param);

    end
    param.delta_2    = 10;
    f            = @(x) ( (1-param.alpha).*2./param.sigma^2.*param.zeta./x./q_0_fun(x,param));
    param.psi    = (1./q_0_fun(param.m_u,param)+integral(f,param.m_l,param.m_u));
    param.mean_x = mean_fun(aux.n,'x',param);
	
    diff=param;
    
end
